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We report on the implementation of Burgers equation as a type-II quantum computation on 
an NMR quantum information processor. Since the flow field evolving under the Burgers equation 
develops sharp features over time, this is a better test of liquid state NMR implementations of type-II 
quantum computers than the previous examples using the diffusion equation. In particular, we show 
that Fourier approximations used in the encoding step are not the dominant error. Small systematic 
errors in the collision operator accumulate and swamp all other errors. We propose, and demonstrate, 
that the accumulation of this error can be avoided to a large extent by replacing the single collision 
operator with a set of operators with random errors and similar fidelities. Experiments have been 
implemented on 16 two-qubit sites for eight successive time steps for the Burgers equation. 

PACS numbers: 



It has been suggested that some classical computa- 
tional problems can be solved by using a hybrid clas- 
sical quantum device, a type II quantum computer 
Such a device is essentially an array of small quantum im- 
formation processors (QIP) sharing information through 
classical channels. NMR has proven to be a useful testbed 
for QIP, and in particular we have shown that a lattice 
of parallel QIPs can be mapped onto a spin system by 
creating a correspondence between the lattice sites and 
spatially distinct spin ensembles. A first proof-of-concept 
for numerically predicting the time-dependent solution 
of classical partial differential equation with dissipative 
terms using our NMR technique was demonstrated for 
the diffusion equation 0, Q . 

One of the most important challenges to implement- 
ing a useful type-II quantum architecture is to avoid the 
accumulation of systematic errors. In the NMR imple- 
mentations to date there are two important sources of 
systematic errors: (1) a linear approximation relating the 
excited magnetization to the Fourier components of the 
shaped RF pulse; and (2) errors from the repeated col- 
lision operators. Here we explore the impact of these 
errors on a simple computation and illustrate a simple 
means of reducing the accumulated error. 

The ensemble nature of the spin system allows us 
to split the sample into a spatial array of lattice sites. 
Well developed methods from magnetic resonance imag- 
ing (MRI) allow us to selectively address the spins in 
each of these sites. Typically the addressing is carried 
out in a space reciprocal to the spatial mapping, called 
fc-space, where k is the wave-number of the correspond- 
ing Fourier components. The fc-space formalism [5| pro- 
vides a recipe for writing a spatially varying spin rotation 
across an ensemble of spins that have been distinguished 
from each other by a magnetic field gradient. The k- 
space formalism is essentially the application of shaped 
radio frequency (RF) pulses in the presence of a linear 
magnetic gradient field as a means of exciting selective 
frequencies. For most studies the full k-space formalism 
is not employed and a linear approximation is invoked. If 



the rotation angle of the shaped pulse is small, then the 
excited magnetization may be accurately calculated only 
to first order in that angle, and the excited magnetiza- 
tion is related to the RF waveform simply by a Fourier 
transform. As a result, the required RF waveform can 
also be determined by taking the inverse Fourier trans- 
form of the desired initial magnetization. This technique 
allows us to encode arbitrary magnetization profiles span- 
ning the various spatial locations in our experiment and 
thereby approximating any desired initial conditions. In 
the previously implemented diffusion equation, higher or- 
der Fourier components of the number density are atten- 
uated by the dynamics and the solution is stable even in 
the presence of substantial accumulated errors. 

To push the development of type-II implementations 
we have chosen to explore the nonlinear Burgers equa- 
tion to test the breakdown for the Fourier approximation. 
Over time, a shock front forms and high spatial frequen- 
cies in the magnetization profile become important and 
it is these high spatial frequencies that we expect to be 
most sensitive to errors. The numerical treatment of the 
QLG algorithm for the Burgers equation |(| therefore of- 
fers a stronger proof of our NMR quantum computing ap- 
proach since the effect of the nonlinear convective term in 
the equation generates a sharp edge as a shock develops 
in time that is not mimicked by spin relaxation, random 
self-diffusion, nor RF inhomogeneities. In addition, we 
demonstrate shock-formation driven by a tunable viscos- 
ity parameter to show that the width of the shock front 
is not determined by implementation imperfections. 

The first-order accurate Fourier approximation was ex- 
pected to be the dominant error source in the NMR 
implementation. However, NMR simulations with con- 
trolled errors shows that the systematic error induced by 
the experimental implementation of the unitary collision 
operator associated with the quantum lattice gas (QLG) 
algorithm is the major challenge. Replacing the single 
collision operator with a set of operators to randomize 
errors allows us to improve the robustness of the imple- 
mentation. 
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A quantum lattice gas is a system of quantum particles 
moving and colliding on a discrete spacetime lattice. This 
quantum particle system is isomorphic to a lattice-based 
qubit system. The mapping is as follows: the probabil- 
ity of a particle residing at a particular lattice node is 
equated to the moduli squared of the probability ampli- 
tude of a qubit at a unique location being in its excited 
quantum state. That is, each spatial location that a par- 
ticle may occupy is mapped onto a qubit associated with 
a unique location. 

The dynamics of evolution in the QLG algorithm can 
be described in three scales, the microscopic, mesoscopic, 
and macroscopic scales. At the microscopic scale, each 
particle has some probability of moving along the lat- 
tice. For example, the particle can move to the right 
or left lattice site in a one-dimentional construction. A 
simplified dynamics allows a particle to change its direc- 
tion of motion (via a collision with another particle) or 
keep moving at constant speed in its original direction 
of motion (streaming). One can describe the dynamical 
behavior of particles at the mesosopic scale by determin- 
ing their occupation probablities on the lattice points; in 
one-dimention, only left-moving and right-moving prob- 
abilities are needed. The mesoscopic dynamical behav- 
ior of the system is modeled by a finite-difference form 
of a quantum Boltzmann equation. Finally, to bridge 
to the macroscopic scale, the occupation probabilities of 
the particles residing at each lattice site are summed to- 
gether to determine the number density. This number 
density quantity defined at each lattice node becomes a 
continuous field at the lattice resolution approaches in- 
finity, which is called the continuum limit. Through a 
Chapman-Enskog purturbation procedure applied in the 
continuum limit, from the quantum Boltzmann equation 
emerges an effective field theory that is parabolic in time 
and space and nonlinear in the number density 0. 

The QLG algorithm is initialized, in the NMR case, by 
encoding the particles' occupation probabilities as a spin- 
magnetization profile. To handle the one-dimensional 
Burgers equation, it is sufficient to use two qubits (two 
spin-| nuclei) per lattice site, where each stores a single 
real valued occupation probability. A room-temperature 
solution of isotopically-labeled chloroform ( lz CHCl^) 
was chosen for implementing the experiments, where the 
hydrogen and the labeled carbon nucleus served as qubits 
f and 2, respectively. The difference of the gyro-magnetic 
ratio of two spins generates widely spaced resonant fre- 
quencies that allows us to address each spin indepen- 
dently. A lattice of QIPs are related to the ensemble 
sample by creating a correspondence between lattice sites 
and spatially dependent positions in the sample. A lin- 
ear magnetic field gradient is used to generate distinct 
spatially-dependent resonant frequencies that we can dis- 
tinguish and modulate by a shaped RF pulse. In this way, 
the magnetic field gradient allows the entire spin ensem- 
ble to be sliced into a lattice of smaller, and individually 
addressable, sub-ensembles. 

The lattice initialization starts by transforming ther- 
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FIG. 1: QLG algorithm implemented in four steps. Three hori- 
zontal lines represent proton spin, carbon spin and fried gradients. 
Both starting magnetizations are encoded in proton channel first 
due to the high signal to noise ratio while decoupled in carbon chan- 
nel to prevent interfering of scalar coupling. The collision operator 
is applied after the initialization. Measurements are also taken in 
two steps in the proton channl followed by data processing in a 
personal computer. 



mal equilibrium states into pseudo-pure states f.8]. The 
equilibrium state is highly mixed and the two nuclear 
spins have unequal magnetizations. Thus, equalization 
of the magnetizations is required prior to creating the 
pseudo-pure state. The dynamical evolution is caused 
by a collision operator (a quantum operation) , and mea- 
surement and streaming (classical operations) according 
to the QLG algorithmic paradigm. The four main sec- 
tions of the NMR implementation of QLG algorithm are 
graphically depicted in Figure ^ 

First, each occupation probablility is mapped onto a 
lattice site as the expectation value of a number opera- 
tor at a spacetime site nAx at time mAt. As a result, 
the initial state of the a th qubit is yj f a (nAx, mAi)\l) + 
i/l - faJnAx, to At) 1 0). The combined the wave func- 
tion for a lattice site is a tensor product over the qubits, 

\^{nAx,mAt)) = y/hh\l\) + VflO- ~ /a)|10> (1) 

+ v / (i-/i)Moi) + V(l-/i)(l-/ 2 )|00) 

In the basis of a two-qubit system, the number opera- 
tors for the occupancy of qubits are defined in terms of 

the singleton qubit number operation n = ^ 







as 



follows: h i = 1 <g> 77 and 77.2 = h (8 1. Therefore, the 
occupation probability is represented as follows: 

f a (?7 Ax, mAt) = {ip(nAx,mAt)\n a \tp(nAx,mAt)). (2) 

The macroscopic scale dynamical quantity of the quan- 
tum lattice gas is the number density, p, defined as the 
sum of the occupancy probablity. The equilibrium occu- 
pation probabilities that we use are 



J a 



eq 



(3) 
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where e a is ±1 for different qubits. 

The initial magnetization is specified by using a RF 
pulse shaped by the Fourier transform of the desired 
magnetization (tranform of the initial number density 
profile). While applying the shaped pulse, a carbon de- 
coupling sequence is performed to prevent the scalar cou- 
pling from interfering with the low power shaped pulses. 
In addition, the ^ pulse, which rotates the information 
from the x-axis to the z-axis, is applied separately just 
after each initialization. This is done to keep the valu- 
able information along the longitudinal direction where 
it will not be affected by the gradient and chemical shift. 
The encoding of initial states on both spins is accom- 
plished in two steps: The initial carbon magnetication is 
recorded on the protons before being transferred to the 
carbons and followed by the initialization of proton mag- 
netication. Furthermore, a short pulse sequence, called 
the clean sequence, is executed after the first swap gate 
to erase the phase distortion that may be caused by the 
decoupling sequence. 

Second, the evolution of f a is governed by the com- 
bined action of the collision operator, measurement and 
streaming. The collision operator is applied to all the 
lattice sites independently, resulting in \ip (nAx)) = 
C\ip(nAx)) 7 for all n. The choice of the particular com- 
ponents of the unitary collision operator determines the 
form of the macroscopic effective field theory (a parabolic 
partial differential equation) and the value of its trans- 
port coefficients (coefficients of the dissipative terms) . A 
general representation of the collision operator for the 
Burgers equation is a block diagonal matrix. This single 
quantum operator is chosen to be 



C = exp 



(a x a y -a y a x ) 



"4.882 y ~ x v v 
which has the following matrix representation 
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The unitary operator C can be decomposed of a se- 
quence of RF pulses and scalar coupling. The product 
operators in the exponent commute with each other, re- 
sulting inC* = exp [-ij^2<j" °f] ex P Hiifes ^ a x ] ■ 
Both terms can be expanded as natural scalar Hamilto- 
nian couplings sandwiched with the appropriate single 
rotations, resulting in 



C = e"^K + ff 
x e *f«+O e 



y , e -i^z °z e - i 4-Wi(' y x -°* ) 
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The exponential terms of single spin rotations are imple- 
mented by 7r/2 and 7r/4 pulses. The exponents of terms 
with <J^ represent the natural internal Hamiltonian 
evolutions with time period 1/2 J, where J is 214Hz. 
Here, the evolution of the internal Hamiltonian is ignored 



while the RF pulse is applied. This approximation leads 
to a systematic error that will accumulate during the 
course of the computation. In general, these errors are 
easy to avoid, but since the purpose of the investigation 
was to explore the sensitivity to accumulated errors we 
did not correct it. The collision operator follows the en- 
coding (Step 2), and it is implemented without magnetic 
field gradients to ensure that all of the sites in the sample 
undergo the same transformation. 

Third, we measure the occupation probabilities. This 
process erases all the superpositions and quantum entan- 
glement that was created by the unitary collision opera- 
tor in the second step. 

The occupation numbers of each spin are ob- 
tained following the collision step by measuring the z- 
magnetization according to the following equation 

f a (n,m) = l - [1+ (^(n,m)\a a z \^(n,m))}. (7) 

Since only a x and a y are observable in our NMR spec- 
trometer, a 7r/2 pulse has been used to bring the z- 
magnetization into the transverse plane. The measure- 
ments are done in two seperate experiments, where a 
SWAP gate is applied to bring the magnetization from 
carbon channel to the proton channel. This SWAP op- 
eration is done because the higher signal-to-noise ratio 
in the proton channel allows us to improve the accuracy 
of our implementation. During the "readout" process 
(Step 3), a week magnetic field gradient is applied to dis- 
tinguish different sites. The observed proton signals are 
digitized and Fourier transformed, allowing us to record 
the spatially-dependent spin magnetization profile. 

Fourth, and last step of the QLG algorithm, we shift 
the f a obtained in the previous step to its nearest neigh- 
bor using a short Matlab program. This step requires 
only classical communication between neighboring sites. 
The time is incremented after this step. Then, we loop 
back to step 1 and update the field of occupation proba- 
bilities over the lattice sites. In this way, we can continue 
to iterate forward in time and make a time-history record 
of the occupation probabilities, which in turn gives us the 
temporal evolution of the number density field. In the 
implementation of the Burger equation, we observed de- 
viations between the numerically predicted data points 
and analytically predicted solutions. These errors can 
be attributed to imperfections in the NMR implementa- 
tion. The major error sources in the NMR implemen- 
tation are known, so to explore the source and relative 
strength of these errors, we have simulated the NMR ex- 
periments. The major error source in this implemen- 
tation is the collision operator, and it is introduced by 
ignoring the scalar coupling between proton and carbon 
during the RF pulses. When applying an RF pulse on 
the proton qubit, the Hamiltonian in the rotating form is 
H = 2nJa^a^' + jnBia x , where B\ is the strength of 
the RF pulse. With the presence of the scalar coupling, 
a small portion of the proton magnetization has been 
transfercd to the carbon qubit. Therefore, the applied 
propagator can be recast as U — U dcsircd U a „ OI . 
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FIG. 2: The growth of the systematic errors due to the collision 
operator in two NMR implementation. The single collision opera- 
tor data (dots) is fit (solid line) with a line of slope 1, which shows 
linear growth of the error. The collision operator data with mod- 
ulated phases (pluses) is the fit with a line of slope 3/4 (dashed 
line). The buildup of the systematic errors has been slowed down 
by proposed method. However, the systematic errors have not been 
totally converted into random errors. 



The error in the collision operator is a systematic er- 
ror that builds up throughout the successive time steps. 
Although this is not the dominant error at the begin- 
ning of the implementation, it eventually dominates the 
first-order error due to the Fourier approximation and 
becomes the dominant issue after just several time step 
interations. Notice that while the reduction of the initial 
magnetization from the Fourier transform is systematic, 
since the magnetization profile is changing the errors are 
not precisely repeated. In the collision operator, how- 
ever, the errors are exactly the same from step to step. 
In addition we expect that the radio frequency inhomo- 
geneity leads to strongly correlated errors in the lattice 
encoding. Hence, we have proposed replacing a single 
collision operator with a set of collision operators that 
have similar fidelity but randomized error terms. 

Since the collision operator for the Burgers equation is 
a zero-order coherence term, the collision operator com- 
mutes with the rotation operator. Therefore, we apply 
a 90° rotation operator to the collision operator at each 
step to mitigate error growth. Consequently, a dramatic 
improvement is observed as shown in Figure 2. On a 
logarithmic plot, the simulation results fit a line with a 
slope of 3/4. If the error terms in the collision operators 
were totally randomized and hence followed a Gaussian 
distribution, the best-fit regression line should have had 
a slope of 1/2. The devitation between our simulation 
data and the ideal Gaussian case indicates residual sys- 
tematic error in the collision operator. In a future study, 
we may use strongly mudulated pulses to randomize the 
error terms. The experimental number densities are over- 
plotted in Figure |21 with the exact analytical solutions. 
Eight successive time steps of the quantum algorithm 



were implemented on 16 two-qubit sites. An improve- 
ment of our p resent experiment a l approach using collision 
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FIG. 3: The experimental data are plotted together with the ana- 
lytical solutions for 8 time steps on a lattice of 16 parallel two-qubit 

QIPs. Viscosity: j . Experimental NMR data (dots) versus 
analytical solution (curves). Randomizing the error terms in the 
collision operator has improved the experimental results dramati- 
cally. 



operators with modulated phases is observed. The agree- 
ment of the data to the analytical solutions is encourag- 
ing and suggests that totally randomizing error terms in 
the collision operator may offer further improvement. 

NMR quantum simulations has provided an alternative 
way to study the NMR spectroscopic implementations. 
From the simulation, we find the major error sources 
are due to imperfect control of the quantum spin system 
and the Fourier approximation associated with setting its 
magnetization profile. Our proposed method for convert- 
ing the systematic errors into random errors is effective. 
The improvement we achieve relative to the previous ex- 
periment is encouraging, and it demonstrates the possi- 
bility of using the same technique in future studies. The 
closeness of the numerical data to the exact analytical 
results for the nonlinear Burgers equation further proves 
the practicality of implementing the QLG algorithm us- 
ing a spatial NMR technique. In addition, although the 
limitation of the Fourier approximation is not dominant, 
the problem of precisely initializing a lattice of QIPs still 
remains an open issue. 
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